ttime=ttime+dt2
pNew=0.0d0
Tnew=0.0d0
Unew=0.0d0
Enew=0.0d0
rhonew=0.0d0

! viscosity 
gradU=0.0d0

do k=1,nInternalFaces

    vtmp=meshSf(1,k)*Usnew(1,k)+meshSf(2,k)*Usnew(2,k)+meshSf(3,k)*Usnew(3,k)
    i=meshFaceOwner(k)
    j=meshFaceNeighbour(k)
    rhonew(i)=rhonew(i)+vtmp*Rhosnew(k)
    stress=0.0d0

    do l=1,3
        do m=1,3
            do m1=1,4
                stress(l,m) = stress(l,m) + gradU(l,m,meshFaces(m1,k)) / 4.0d0 ! kg/m3 * m2/s2 = kg / (m * s2)
            end do
        end do
    end do

    DivVStress = 0.0d0

    do l = 1,3
        do m = 1,3
            DivVStress = DivVStress + meshSf(l,k) * Usnew(m,k) * stress(l,m)
        end do
    end do

    Enew(i)=Enew(i)+vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) + DivVStress
    Unew(1:3,i)=Unew(1:3,i)+vtmp*Rhosnew(k)*Usnew(1:3,k)+(Psnew(k)-1d5)*meshSf(1:3,k) +meshSf(1,k)*stress(1,1:3) &
                           +meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)
    rhonew(j)=rhonew(j)-vtmp*Rhosnew(k)

    Enew(j)=Enew(j)-vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) - DivVStress
    Unew(1:3,j)=Unew(1:3,j)-vtmp*Rhosnew(k)*Usnew(1:3,k)-(Psnew(k)-1d5)*meshSf(1:3,k) -meshSf(1,k)*stress(1,1:3) &
                           -meshSf(2,k)*stress(2,1:3)-meshSf(3,k)*stress(3,1:3)
    ! [Unew] = m2 * stress = m2 * kg / (m * s2) = kg * m / s2

end do		

! cycle on boundary faces

do k=nInternalFaces+1,nFaces

    vtmp=meshSf(1,k)*Usnew(1,k)+meshSf(2,k)*Usnew(2,k)+meshSf(3,k)*Usnew(3,k)
    i=meshFaceOwner(k)
    stress=0.0d0

    do l=1,3
        do m=1,3
            do m1=1,4
                stress(l,m)=stress(l,m) + gradU(l,m,meshFaces(m1,k))/4.0d0
            end do
        end do
    end do

    DivVStress = 0.0d0
 
    do l = 1,3
        do m = 1,3
            DivVStress = DivVStress + meshSf(l,k) * Usnew(m,k) * stress(l,m)
        end do
    end do

    rhonew(i)=rhonew(i)+vtmp*Rhosnew(k)
    Enew(i)=Enew(i)+vtmp*(Rhosnew(k)*Esnew(k)+Psnew(k)) + DivVStress
    Unew(1:3,i)=Unew(1:3,i)+vtmp*Rhosnew(k)*Usnew(1:3,k)+ (Psnew(k)-1d5)*meshSf(1:3,k) +meshSf(1,k)*stress(1,1:3) &
                +meshSf(2,k)*stress(2,1:3)+meshSf(3,k)*stress(3,1:3)

end do

do i=1,nCells

    rhonew(i)=rho(i)-rhonew(i)/meshV(i)*dt2
    Enew(i) = (E(i)*rho(i)-Enew(i)/meshV(i)*dt2)/rhonew(i) ! E(i) !

    if( SourcePresent ) then
        if ( i .eq. iCellSource ) then
            call CalcSource(ttime,Q)
            Q = Q*dt2  / meshV(i)
            Enew(i) = Enew(i) + Q / rhonew(i)
        end if
    end if

    Unew(1:3,i)=(U(1:3,i)*rho(i)-Unew(1:3,i)/meshV(i)*dt2)/rhonew(i)
    ! [Unew] = kg * m / s2 / m3 * s / (kg/m3) = 1 / s2 / m2 * s * m3 = m / s

    pnew(i)=(gam-1.0d0)*(Enew(i)-(Unew(1,i)**2+Unew(2,i)**2+Unew(3,i)**2)/2.0d0)*rhonew(i)
    Tnew(i)=pnew(i)/(rhonew(i)*cv*(Gam-1.0d0))

    ! determaine Currant's Time Step
    vtmp=DbCells(i)/(dsqrt(Unew(1,i)**2+Unew(2,i)**2+Unew(3,i)**2)+dsqrt(gam*pnew(i)/rhonew(i)))
    if (vtmp.lt.dtmin) dtmin=vtmp
end do

